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Evolution of Stellar Collision Products in Globular Clusters — II. OfF-axis 

Collisions. 

Alison Sills^, Joshua A. Faber^, James C. Lombardi, Jr.^, Frederic A. Rasio^, Aaron R. Warren^ 

ABSTRACT 

We continue our exploration of collisionally merged stars in the blue straggler re- 
gion of the color-magnitude diagram. We report the results of new SPH calculations 
of parabolic collisions between two main-sequence stars, with the initial structure and 
composition profiles of the parent stars having been determined from stellar evolution 
calculations. Parallelization of the SPH code has permitted much higher numerical 
resolution of the hydrodynamics. We also present evolutionary tracks for the resulting 
collision products, which emerge as rapidly rotating blue stragglers. The rotating col- 
lision products are brighter, bluer and remain on the main sequence longer than their 
non-rotating counterparts. In addition, they retain their rapid rotation rates throughout 
their main sequence lifetime. Rotationally-induced mixing strongly affects the evolution 
of the collision products, although it is not sufficient to mix the entire star. We discuss 
the implications of these results for studies of blue straggler populations in clusters. 
This work shows that off-axis collision products cannot become blue stragglers unless 
they lose a large fraction of their initial angular momentum. The mechanism for this 
loss is not apparent, although some possibilities are discussed. 

Subject headings: blue stragglers ~ globular clusters: general - stars: evolution - hy- 
drodynamics - stars: rotation 

1. INTRODUCTION 

Globular clusters have, until recently, been studied under the auspices of either stellar evolution 
or stellar dynamics. Stellar evolutionists often treat clusters as collections of isolated stars all with 
the same age and primordial chemical abundances, while stellar dynamicists often treat clusters 

^Department of Astronomy, The Ohio State University, 140 W. 18th Ave., Columbus, OH 43210; 
asills@astronomy.ohio-state.edu 

^Department of Physics, MIT, 77 Massachusetts Ave, Cambridge, MA 02139; jfaber@mit.edu, rasio@mit.edu 

■^Department of Physics and Astronomy, Vassar College, 124 Raymond Ave., Box 562, Poughkeepsie, NY 12603; 
lombardi@vassar.edu, aawarrenOO@alum.vassar.edu 



- 2 - 



as a system of point masses moving under the influence of their mutual gravities. However, it is 

becoming clear that the dynamics of a cluster can affect its populations, and vice versa (Bailyn 
1995; Portegies Zwart et al. 1999; Joshi, Nave &: Rasio 2000). Therefore, in order to model globular 
clusters effectively, we need to consider both stellar evolution and stellar dynamics, and especially 
their influence on each other. 

One situation of particular interest in which stellar dynamics and stellar evolution interact 
involves direct stellar collisions. These collisions, especially those involving binary stars, can have a 
significant impact on the energy budget of a cluster. In particular, they can destroy and change the 
properties of hard primordial binaries. Dynamical interactions involving hard primordial binaries 
are thought to be the most important source of heating that can support a cluster against core 
collapse (see, e.g., Goodman & Hut 1989; Gao et al. 1991; Rasio 2000). In addition, it has been 
suggested that, through the indirect heating by mass loss associated with the stellar evolution of 
collision products in the cluster core, collisions could provide a mechanism for supporting a cluster 
against gravothermal collapse after exhaustion of its primordial binaries (Goodman &; Hernquist 
1991). The products of stellar collisions between main-sequence stars have the additional useful 
property that they are potentially observable in the color-magnitude diagrams of clusters as blue 
stragglers. Since the population of blue stragglers is determined by the dynamics of the cluster, we 
can use the blue stragglers to investigate the dynamical state of the cluster. For example, Ferraro 
et al. (1999) have suggested that the large population of blue stragglers in M80 implies that the 
cluster is close to core collapse. In order to use blue stragglers to study globular clusters, however, 
we need to understand clearly the evolution and observable properties of collision products. 

In previous work (Sills et al. 1997, hereafter Paper I), we used the results of Smoothed Par- 
ticle Hydrodynamics (SPH) simulations of stellar collisions as starting models for stellar evolution 
calculations, and produced detailed evolutionary models for head-on (i.e., zero impact parameter) 
collision products in globular clusters. However, collisions between stars are, of course, never exactly 
head-on. Lombardi, Rasio &: Shapiro (1996) calculated collisions with non-zero impact parameters, 
and some of the resulting off-axis collision products were rotating close to their break-up speeds 
when they returned to dynamical equilibrium. The effects of rotation on the subsequent evolution 
of these collision products will be substantial. These effects have never before been included in any 
theoretical study of blue straggler evolution. 

In this paper, we present four new SPH simulations of off-axis collisions, and use the results to 
initiate stellar evolution calculations of the rotating collision products. In §2, we outline the method 
and results of the SPH calculations. We summarize our technique for combining the results of SPH 
simulations of stellar collisions with detailed evolutionary models, and discuss the modifications we 
have made for rotating stars. §3 presents the implementation of rotation used in the Yale Rotational 
Evolution Code (YREC), and the results of the evolution calculations. We discuss the implications 
of our results for blue straggler creation mechanisms and for using blue stragglers as tracers of the 
dynamical evolution of globular clusters in §4. 
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2. HYDRODYNAMIC SIMULATIONS OF OFF-AXIS STELLAR COLLISIONS 

2.1. Implementation of the SPH Code 

Our hydro dynamic simulations were performed using a new, parallel version of the SPH code 
described in Paper I. This code was originally developed by Rasio (1991) specifically for the study of 
hydrodynamic stellar interactions such as collisions and binary mergers (see, e.g., Rasio &: Shapiro 
1991, 1992, 1994). SPH is a Lagrangian hydrodynamics technique in which pressure gradient forces 
are calculated by kernel estimation directly from particle positions (for a review, see Monaghan 
1992). In SPH simulations we solve the equations of motion of a large number N of Lagrangian 
fluid particles moving under the influence of hydrodynamic forces as well as the fluid self-gravity 
(in this paper A'' is as large as 1.05 x 10^). For a basic overview of the SPH equations, see Rasio & 
Lombardi (1999). 

Local densities and hydrodynamic forces at each particle position are calculated by smoothing 
over A^jv nearest neighbors. The size of each particle's smoothing kernel is evolved in time to keep 
Nn close to a predetermined optimal value. For the high-resolution calculations in this paper the 
optimal number of neighbors was set at = 100. Neighbor lists for each particle are recomputed 
at every iteration using a linked-list, grid-based parallel algorithm. Note that neighborhood is not 
a symmetric property, as one particle can be included in another's neighbor list but not vice versa. 
Pressure and artificial viscosity (AV) forces are calculated by a "gather-scatter" method, looping 
over all particles in turn, gathering the force contribution on each particle from its neighbors, and 
scattering out the equal and opposite force contribution on the neighbors from each particle. 

The fluid self-gravity in our code is calculated by an FFT-based convolution method. The 
density field is placed on a 3D grid by a cloud-in-cell method, and convolved with a kernel function 
calculated once during the initialization of the simulation. We use zero-padding of our grids to 
obtain correct boundary conditions for an isolated system, at the expense of memory storage. 
Gravitational forces are calculated from the gravitational potential by finite differencing on the 
grid, and then interpolated for each particle using the same cloud-in-cell assignment. 

A number of variables are associated with each particle i, including its mass rrtj, position rj, 

velocity v,-, cntropic variable Ai and numerical smoothing length hi. The entropic variable ^ is a 
measure of the fluid's compressibility and is closely related (but not equal) to specific entropy: for 
example, both A and specific entropy are conserved in the absence of shocks. For convenience, we 
will refer to A as entropy. In this paper we adopt an equation of state appropriate for a monatomic 
ideal gas: that is, the adiabatic index 7 = 5/3 so that Pi = Aip^J^ , where Pi and pi are the density 
and pressure of particle i. 

The rate of increase, due to shocks, for the entropy Ai of particle i is given by 




(1) 
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where the summation is over neighbors and Wij is a symmetrized smoothing kernel. We adopt the 
form of AV proposed by Balsara (1995): 

where 

/^-CTF^-pT^^^ it(v,-v,).(r,-r,)<0_ 
if (v,-v,)-(r,-r,)>0 

The terms hij and Cij are, respectively, the average smoothing length and sound speed associated 
with particles i and j. Here fi is the so-called form function for particle i, defined by 

f. = l^-"^!' (4) 

\V -vli + lV x-v\i + rfci/hi' ^ ^ 

where 

(V ■ v), = i ^ rrij {Vj - V,) • ViWij (5) 
^- 

and ^ 

(V X v)i = — ^ mj(vi - Vj) X VjWy. (6) 

j 

The function /j acts as a switch, approaching unity in regions of strong compression (|V • v|j >> 
|V X v|i) and vanishing in regions of large vorticity (|V x v|j >> |V • v|i). Consequently, this AV 
has the advantage that it is suppressed in shear layers. In this paper we use 77^ = 0.01, rj' = 10~^ 
and a = /3 = 7/2. This choice of AV treats shocks well, while introducing only relatively small 
amounts of numerical viscosity (Lombardi et al. 1999). 

The code was parallelized using MPI (the message passing interface), to run efficiently on 
multiple processors. Simpler methods, such as those in HPF, prove impossible to use for SPH 
loops, since the "gather-scatter" technique, which enforces Newton's Second Law, disallows the use 
of DO INDEPENDENT loops. In our code, the lists of neighbors are split and divided up evenly 
between the processors. All loops over SPH particles containing sums over quantities involving 
neighbors, including the AV, are distributed among the processors, summed in part, and collected 
at the end of the relevant subroutine by an MPI_ALLGATHER command. The gravitational force 
calculations have been parallelized using the rfftwnd-mpi package of the FFTW library (Frigo & 
Johnson 1997). All 3D grids used by gravity routines are split up and distributed evenly among 
processors in the z-direction, which saves a considerable amount of computer memory. Force 
contributions from grid cells contained within each process are calculated in parallel, and gathered 
at the end of the subroutine by an MPLALLGATHER command. In benchmarking tests, we 
have found that our parallel code scales very well when using up to 32 processors on a distributed 
shared-memory supercomputer. For example, doubling the number of processors up to 32 results in 
an increase of total CPU time by no more than 10% on the SGI/Cray Origin2000 supercomputer. 
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Both the hydrodynamics and gravitational subroutines show similar speed-up when increasing the 
number of processors. 

Our simulations employ equal mass SPH particles, both to keep the resolution high in the 
stellar cores and to minimize spurious mixing during the simulation; see Lombardi et al. (1999) 
for a discussion of spurious transport induced by unequal mass particles. Due to extreme central 
densities in the parent stars, the Courant stability condition requires exceedingly small timesteps. 
Consequently, an elapsed physical time of one hour in our high resolution calculations requires 
roughly 1000 iterations and 100 CPU hours on the SGI/Cray Origin2000 supercomputer at NCSA. 
The three high resolution calculations presented in this paper each required a total of ~ 1000 CPU 
hours. 



2.2. Initial Conditions 

All of our hydrodynamic simulations use parent star models whose characteristics are based 
on the results of YREC calculations, as discussed in Sills & Lombardi (1997). In particular, we 
evolved (non-rotating) main-sequence stars of total mass M = 0.6 and Q.SMq with a primordial 
helium abundance Y = 0.25 and metallicity Z = 0.001 for 15 Gyr, the amount of time needed to 
exhaust the hydrogen in the center of the O.SMq star. The thermodynamic and structural profiles 
of the parent stars are shown in Figure 1. The total helium mass fractions for the 0.6 and O.SMq 
parent stars are 0.286 and 0.395, and their radii are 0.52i?Q and O.96i?0, respectively. From the 
pressure and density profiles of these models, we compute the entropy profile and assign values of A 
to SPH particles accordingly. In addition, the chemical abundance profiles of 15 different elements 
are used to set the composition of the SPH particles. To minimize numerical noise, each parent 
star's SPH model is relaxed to equilibrium using an artificial drag force on the particles, and then 
these relaxed models are used to initiate the collision calculations. 

The stars are initially non-rotating and separated by 5 Rto, where Rto = O.96i?0 is the 
radius of a turnoff star. The initial velocities are calculated by approximating the stars as point 
masses on an orbit with zero orbital energy (i.e., we assume the parabolic limit, appropriate for 
globular clusters) and a pericenter separation r^. A Cartesian coordinate system is chosen such 
that these hypothetical point masses of mass Mi and M2 would reach pericenter at positions 
Xi = (— 1)*(1 — Mi/{Mi + M2))rp, Vi = Zi = 0, where i = 1,2 and i = 1 refers to the more massive 
star. The orbital plane is chosen to be z = 0. With these choices, the center of mass resides at the 
the origin. 

2.3. Hydrodynamic Results and Conversion to YREC Format 

Table 1 summarizes the parameters and results of all of our SPH simulations. The first column 
gives the name by which the calculation is referred to in this paper; we use lower case letters to 
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distinguish the present calculations from the corresponding calculations involving polytropic parent 
stars presented in Lombard!, Rasio & Shapiro (1996). The second and third columns give the masses 
Ml and M2 of the parent stars. Column (4) gives the ratio rp/{Ri +R2), where rp is the pcriccnter 
separation for the initial orbit and Ri + R2 is the sum of the two (unperturbed) stellar radii. This 
ratio has the value for a head-on collision, and ~ 1 for a grazing encounter. Column (5) gives 
the number of SPH particles. Column (6) gives the final time at which the calculation was 
terminated. Column (7) gives the total angular momentum of the merger remnant. Column (8) 
gives the ratio T/|VF| of rotational kinetic energy to gravitational binding energy of the (bound) 
merger remnant in its ccntcr-of-mass frame at time tf. Column (9) gives the mass of the merger 
remnant (the gravitationally bound fluid). Columns (10) and (11) give the velocity components Vx 
and Vy for the merger remnant's center of mass at time tf in the system's center-of-mass frame. 
Since the amount of mass ejected during a parabolic collision is small, the merger remnant never 
acquires a large recoil velocity. The case e, f and k simulations implemented Njsf = 100 neighbors 
per particle, while case j' implemented A^jv = 32 neighbors. Case j' is a low resolution calculation 
of a nearly head-on collision. 

Figure 2 illustrates the dynamical evolution for turnoff main-sequence star (Mi = 

O.8M0) collides with a slightly less massive star (M2 = 0.6Mq). The parabolic trajectory has a 
pericenter separation Vp = 0.25(i?i +-^2)- The first collision at time t ~ 2 hours disrupts the outer 
layers of the two stars, but leaves their inner cores essentially undisturbed. The two components 
withdraw to apocenter at t ~ 3.2 hours before colliding for the second, and final, time. The merger 
remnant is rapidly and differentially rotating, and shear makes the bound fiuid quickly approach 
axisymmetry within a few dynamical timescales. 

Figure 3 shows contours of the density p, the entropy A, and the z— component of the specific 
angular momentum j = ojr'^yi for each collision product in a slice containing its rotation axis. 
Here r^yi is the cylindrical radius measured from the rotation axis. The contours shown represent 
an average over the azimuthal angle, as the product approaches axisymmetry at the end of the 
simulation. By definition, the collision products are not barotropes since the specific angular 
momentum j (and hence u) has a distinct dependence on z. For a rotating star, stable thermal 
equilibrium requires du/dz = in chemically homogeneous regions, as discussed in §3.2. From the 
set of specific angular momentum contours presented in Figure 3, it is therefore evident that the 
outer layers of the collision products (where the composition gradient is small) are not in thermal 
equilibrium. The shape of the specific angular momentum contours is due to the tendency of 
dynamical shear to equalize the angular velocity uj on surfaces of constant density (Pinsonneault 
et al. 1989), so that j increases with z at fixed Vcyi- 

To convert our three-dimensional hydrodynamic results into the one-dimensional YREC for- 
mat, we begin by averaging the entropy and specific angular momentum values in ~ 25 bins in 
enclosed mass fraction. The resulting profiles are extrapolated over the outermost ~5% of the rem- 
nant mass, since in practice the simulations must be terminated before all bound fluid has fallen 
back to the remnant's surface (for a discussion of how the values of the final time tf are chosen, see 
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§2.3 of Lombardi, Rasio & Shapiro 1996). The specific angular momentum profile is then renor- 
malized (by a factor of close to unity), to ensure that the extrapolation does not modify the total 
angular momentum of the remnant. With the A and j profiles given, the structure of the rotating 
remnant is uniquely determined in the formalism of Endal & Sofia (1976), by integrating the gen- 
eral form of the equation of hydrostatic equilibrium [see their eq. (9)]. To do so, we implement an 
iterative procedure in which initial guesses at the central pressure and angular velocity are refined 
until a self-consistent YREC model is converged upon. The temperature profile is calculated from 
the ideal gas equation of state, while the luminosity profile is calculated according to the equation 
of radiative transport. This first approximation for the structure of the star was used as a starting 
model for the stellar evolution code. 

Although the pressure, density, and angular velocity profiles of the resulting YREC model 
cannot also be simultaneously constrained to equal exactly the corresponding (averaged) SPH pro- 
files, the differences are slight (see, for example, the uj profiles in Figure 4). Indeed, the subsequent 
evolutionary tracks are not significantly different for remnant models that use the form of the to 
profile taken directly from SPH. The decision to constrain the j profile, as opposed to the angular 
velocity uj profile, guarantees that the total angular momentum of the initial YREC model agrees 
with that from the SPH results. Furthermore, the SPH uj profile can be slightly less accurate than 
the J profile, since the latter is less sensitive to the resolution of the calculation. 



3. EVOLUTION OF THE ROTATING COLLISION PRODUCTS 

3.1. Implementation of the Stellar Evolution Code 

Our stellar evolution calculations are performed with YREC. YREC is a one-dimensional 
code, in which the star is divided into shells along surfaces of constant gravitational plus rotational 
potential. The code solves the equations of stellar structure with the Henyey technique, and 
follows the rotational evolution with the formalism of Endal & Sofia (1976). Guenther et al. (1992) 
give a detailed description of the physics implemented in the evolution code. We used the same 
opacities, equation of state and model atmospheres as in Paper I. The free parameters in the code 
(the mixing length and parameters that set the efficiency of angular momentum transport and 
rotational chemical mixing) are set by calibrating a solar mass and solar metallicity model to the 
sun. 

Rotation is treated by evaluating physical quantities on equipotential surfaces rather than 
the spherical surfaces usually used in stellar models. The hydrostatic equilibrium and radiative 
transport equations contain terms that account for the lack of spherical symmetry in the rotating 
star. A number of rotational instabilities that transport angular momentum and material within 
the star are followed, including dynamical shear (Pinsonneault et al. 1989), meridional circulation 
(von Zeipel 1924), secular shear (Zahn 1974), and the Goldreich-Schubert-Frickc (GSF) instabil- 
ity (Goldreich & Schubert 1967; Pricke 1968). Angular momentum transport and the associated 
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chemical mixing are treated as diffusion processes, witfi diffusion coefficients tfiat account for eacfi 

active mechanism within unstable regions of the star. The diffusion coefficients arc proportional to 
circulation velocities, which have been estimated by Endal & Sofia (1987). In addition to the inter- 
nal rearrangement of angular momentum, angular momentum also can be drained from an outer 
convection zone through a magnetic wind, using the formalism given by Chaboyer, Demarque k, 
Pinsonneault (1995). The Endal &; Sofia scheme is valid across a wide range of rotation rates, for a 
restricted class of angular momentum distributions. This scheme requires that the angular velocity 
is constant on equipotential surfaces, which does not allow for modeling of latitude-dependent an- 
gular velocity profiles. See Meynet &: Maeder (1997) for a detailed discussion of the validity of this 
approach. For a detailed description of the implementation of rotation in YREC, see Pinsonneault 
et al. (1989). 



3.2. Evolutionary Results 

Although in dynamical equilibrium, the nascent collision products are initially far from thermal 
equilibrium and contract on a thermal timescale. During this stage, their luminosity is predomi- 
nately provided by the conversion of gravitational potential energy, analogous to the situation in 
pre-main-sequence stars. After the contraction, central hydrogen burning quickly becomes signifi- 
cant, and we refer to this point as the zero age main sequence (ZAMS). Note, however, that unlike 
normal ZAMS stars, these stars are far from being chemically homogeneous. Many qualitative fea- 
tures of the subsequent evolution are similar to that of normal intermediate mass stars; for example, 
it is possible to define the subgiant and giant stages of evolution. Nevertheless, the timescales for 
this evolution, and even the actual evolutionary tracks, can be substantially different from those of 
normal stars, because of the initial structure of the collision product and because of the effects of 
the high rotation rate. 

The initial total angular momentum of the collision products (Table 1) is as much as 10 
times larger than that of normal stars of comparable mass at the pre-main-sequence birth line. 
The stellar collision products differ substantially from normal pre-main-sequence stars in another 
significant way. Normal pre-main-sequence stars have large convective envelopes, and therefore can 
lose substantial amounts of angular momentum via a magnetic wind. The sun has lost about 99% of 
its initial angular momentum, for example (Pinsonneault et al. 1989). However, the stellar collision 
products are not convective at the end of the SPH simulation, nor do they develop convection zones 
during the thermal readjustment and main-sequence stages. Therefore, the collision products do 
not lose any angular momentum by a magnetic wind until the subgiant branch. 

Rotational instabilities transport angular momentum from the center of the star to the sur- 
face approximately on an Eddington-Sweet timescale (the thermal timescale divided by the ratio 
of centrifugal acceleration to gravitational acceleration). For the rapidly rotating stellar collision 
products, this timescale is on the order of only 10^ years for the stars as a whole, and is even 
shorter in the outer layers. The dominant mixing mechanism in these stars is the GSF instability. 
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which addresses axisymmetric perturbations: A structure in which the angular velocity uj is con- 
stant on cylinders is the only thermally stable configuration in regions of homogeneous chemical 
composition. Since dynamical shear acts to remove any latitude dependence of the angular veloc- 
ity (see Pinsonneault et al. 1989), regions that are subject to the GSF instability evolve towards 
rigid rotation. In the case of our stellar collision products, the angular velocity profile is initially 
a steeply decreasing function of enclosed mass fraction (Figure 4). Therefore, the GSF instability 
serves to raise the surface rotation rate and decrease the interior rotation rate, flattening the co 
profile. Since the timescale for this to occur is shortest at the surface, the uj profile begins to fiatten 
there, and then the fiattened profile penetrates inwards. The star has a very large reservoir of 
angular momentum, and this redistribution of angular momentum can cause the surface of the star 
to rotate faster than its breakup velocity. 

An obvious solution here would be simply to remove the mass from the calculation as it becomes 
unbound to the star. Unfortunately, under this scheme, rapidly rotating stellar collision products 
could be completely disrupted. The total angular momentum of the star is so large that if the star 
were rotating as a solid body, its rotational velocity would be more than its breakup velocity. This 
is true immediately following the collision, when the star has its largest radius, and is even more 
obvious as the star collapses to the main sequence. The breakup velocity is proportional to 
so it increases as the radius R decreases. However, the surface rotation rate is proportional to R~^, 
and increases faster than the breakup velocity as R decreases. As the GSF instability works to 
fiatten the rotation profile further and further into the center of the star, more of the star becomes 
unbound, and eventually the entire star can be spun apart. This poses a problem for the theory 
that blue stragglers are created via stellar collisions: either blue stragglers are not created through 
physical off-axis stellar collisions, or some mechanism(s) can remove angular momentum from the 
star on short timescalcs. Possible mechanisms will be discussed in §4. 

For now wc will simply postulate that the angular momentum of the collision products can be 
reduced. In order to study the effects of rotation on stellar collision products, we performed some 
experiments on case k {O.GMq + 0.6MQ,rp = 0.25(i?i -|- -^2))- We artificially decreased the initial 
rotation velocity of each shell by a large factor before evolving the model. In particular, we divided 
the rotation velocities by 1000, 100, 10, 5 and 2, which reduced their total angular momentum 
by almost the same factors. In this last case, the star still had enough angular momentum that 
it began to rotate faster than its breakup velocity. The evolutionary tracks for the other angular 
momentum values are shown in Figure 5. The tracks with the least angular momentum are not 
noticeably different from the evolutionary track of the corresponding non-rotating model. As the 
total angular momentum of the star increases, the star spends more time on the main sequence due 
to two effects. First, rotation provides an additional source of support for the star, so that it does 
not have to produce as much pressure through nuclear burning to prevent its collapse. Second, 
rotationally induced mixing moves hydrogen into the core of the star, increasing the amount of 
available fuel. 

We also evolved the product of the low resolution SPH simulation of a barely off-axis collision 
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(case j'). In this situation, two 0.6 Mq parent stars collide with an impact parameter of only 
Tp = 0.02i?, where R = 0.52i?Q. The SPH results are used directly as the initial conditions for the 
evolution code as described in §2.3, without rescaling the angular velocity profile. This collision has 
sufficiently little angular momentum that even if the collision product rotated as a solid body, it 
would rotate slower than its breakup velocity. We used this case to confirm that the results of SPH 
simulations could indeed be evolved directly, and to see what the effect of even a barely off-axis 
collision was on the subsequent evolution of the collision product. The results are shown in Figure 
6. For comparison, a non-rotating version of the same star is also shown. The non-rotating star 
was evolved using the same starting model as the rotating case, but with zero angular momentum. 
The effects of rotation, even for a collision that is barely off-axis, are remarkable. This collision 
product has an initial total angular momentum of 2.1 x 10^° g cm^ s~^, comparable to that of a pre- 
main-sequence star of the same mass. Since the collision product does not lose angular momentum 
during its evolution through the main-sequence stage, its rotational velocity remains quite large 
(272 km s^^ at the turnoff). Substantial rotational mixing occurs, prolonging its main-sequence 
lifetime to 3.7 Gyr (compared to 1.8 Gyr for the equivalent non-rotating star). 

The evolutionary tracks for the case e, f and k collision products are shown in Figure 7. Here 
we have reduced the initial angular velocities by a factor of 5, in order to prevent the stars from 
rotating faster than breakup. None of the other profiles are changed initially, although YREC 
converges on a solution which satisfies the rotational equations of stellar structure. Some details of 
the various evolutionary tracks are given in Table 2. Cases e and f involve the same parent stars 
(0.6 Mq and 0.8 Mq) but have different pericenter separations, namely 0.25 and 0.5 (i?i -|- R2), 
respectively. Since case f has significantly more angular momentum than case e, much more mixing 
occurs, and the evolutionary tracks are substantially different. The track of the case f collision 
product also shows the most striking differences from the track of a normal, slowly rotating star. 
Although not even the case f collision product is fully mixed, hydrogen is carried into the core 
nearly as fast as it is depleted by nuclear burning, extending the main-sequence lifetime in case f 
by a factor of 5 over that of case e. 

In addition to transporting angular momentum and carrying hydrogen into the core, rotational 
instabilities also carry helium to the surface, affecting the surface opacity and causing the star to 
become continuously brighter and bluer while on the main sequence. The most dramatic example 
of this effect can be seem in the case f track of Figure 7, which reaches a luminosity of 300Lq and a 
surface temperature of 2.5 x 10"^ K while still on the main sequence. In this case, the helium mass 
fraction at the surface reaches an incredible Y = 0.85. Figure 8 shows the surface abundance of 
helium as a function of time for the collision products featured in this paper. 

The differences between the case e and case f tracks, which correspond to collision products 
of nearly the same mass, help to emphasize that the total mass of a collision product is not the 
primary factor in determining its relative position on a Hertzsprung-Russell (HR) diagram, unlike 
the situation for normal main-sequence stars. As a second example, consider the collision product 
of case k, which has a mass of only I.I6M0, but which is nevertheless brighter and bluer at turnoff 
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than the 1.36M0 product of case e. As discussed in Paper I and Sills Sz Lombardi (1997), a collision 
that involves a turnoff mass star (such as in cases e and f) will produce a blue straggler with a 
hydrogen-depleted core, since the core of the turnoff star sinks to the center of the collision product. 
Therefore, the main-sequence lifetime for the higher mass products is generally shorter. In lower 
mass coUisional blue stragglers, there is then more time for helium to be mixed to the surface, and 
the surface helium abundance can grow larger than in the higher mass collision products (compare 
cases e and k in Figure 8). As a result, a lower mass collision product can become brighter and 
bluer than one of higher mass. Consequently, blue straggler masses that have been determined 
from color-magnitude diagram position could be in error if the blue stragglers are rapidly rotating. 

Figure 9 shows the surface rotational velocity as a function of time for the four cases featured 
in this paper. The ages of the ZAMS, the turnoff, and the base of the giant branch are marked 
in each of the panels. All four cases show similar behavior: the rotational velocity adjusts rapidly 
as the star contracts to the main sequence. Once on the main sequence, the surface rotational 
velocity remains approximately constant. As the star becomes a subgiant, it develops a surface 
convection zone, and the standard angular momentum loss mechanism (via a magnetic wind) begins 
to operate. The loss rate is proportional to the surface angular velocity to the third power, so these 
rapidly rotating stars very efficiently shed most of their angular momentum before reaching the 
giant branch. 

4. SUMMARY AND DISCUSSION 

We have modeled blue stragglers created in off-axis collisions between cluster main-sequence 
stars. Since even a minutely off-axis collision can produce a rapidly rotating star, most collision 
products cannot remain gravitationally bound as they contract to the main sequence unless they 
somehow lose angular momentum. The standard angular momentum loss mechanism (via a mag- 
netic wind) docs not operate in collisional blue stragglers, since magnetic braking is effective only 
in stars with surface convection zones deep enough to support a dynamo. The stellar collision prod- 
ucts lack a surface convection zone, even during their collapse to the main sequence. Therefore, 
either blue stragglers cannot be created by stellar collisions, or some other mechanism for removing 
angular momentum must occur. 

Angular momentum can be lost from stars in a variety of ways. Stars more massive than these 
collision products have substantial mass loss rates during their main-sequence evolution, which can 
be on the order of 10~^Mq yT~^ for stars of 3 Mq (Lamers 1981). Such a high mass loss rate 
will remove a substantial amount of angular momentum from a star. The mass loss in these stars 
takes the form of a line radiation-driven wind (Castor, Abbott & Klein 1975), which requires high 
luminosities in order to remove a large amount of mass. This mechanism may be of importance in 
extremely bright collision products (such as in case f), although low luminosity blue stragglers are 
unlikely to be significantly affected. 
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If the collision product is a member of a binary system, it can lose angular momentum through 
the orbital evolution of the system. Angular momentum loss is usually discussed in the context of 
contact or near-contact binary systems, and requires that one of the stars has a stellar wind and 
a co-rotating magnetic field (Stcpien 1995). It is possible to have a stellar collision in which the 
product has a binary companion (through an encounter of a single star with a binary system (Davies 
1995), or two binary systems with each other (Bacon, Sigurdsson, &; Davies 1996). However, it is 
not clear that the resulting binary system will be in contact, and it is highly unlikely that all blue 
stragglers formed in such a fashion. 

The most plausible mechanism for removing angular momentum from these stellar collision 
products is disk-locking. In low mass pre-main-sequence stars, the torque placed on a star through 
a magnetic field coupled to a disk of material at a few stellar radii is enough to slow the rotation 
of these stars by as much as an order of magnitude (Sills, Pinsonneault k. Terndrup 2000). This 
mechanism requires the presence of both a magnetic field and a disk of material around the star. 
The SPH simulations of the stellar collisions show that typically a few percent of the mass of the 
system is not bound to the final collision product. However, this mass will leave the system on a 
dynamical time (on the order of hours) and cannot become the disk. A more reasonable source of 
disk mass is the outer layers of the star that are spun off during the initial contraction on a thermal 
timescale. Since the amount of mass necessary for such a disk is on the order of a few hundredths 
of a solar mass, the disk-locking could begin very soon after the collision, and remove the angular 
momentum quite quickly. A simple calculation shows that this scenario is not unreasonable. We 
assume that angular momentum is carried away by mass which is coupled to the surface rotation 
rate of the star down to the Alfven radius (Konigl 1991). In order to lose 80% of the star's angular 
momentum (in case f), we need to lose only ~ O.IM©. The other ingredient in this scenario is the 
presence of a magnetic field. We know that low mass stars have substantial magnetic fields. It is 
not clear what happens to those fields after the stellar collision, and whether they will be coherent 
enough to lock the star to a disk. More work is necessary to determine if this or any other angular 
momentum loss mechanism occurs in stellar collision products. 

We assumed that the stellar collision products would lose a large fraction (up to ~80%) of 
their initial angular momentum early in their evolution. We then continued their evolution through 
the main sequence to the base of the giant branch. The rotating stars are bluer and brighter than 
their non-rotating counterparts, and have longer main-sequence lifetimes. Both of these changes 
to the evolutionary tracks will have an effect on the derived blue straggler distributions in the 
color-magnitude diagram, and therefore will change our conclusions about the dynamical state of 
globular clusters. 

Since rotating blue stragglers live on the main sequence longer, we expect that their distribu- 
tions will be more populated near the main sequence than those presented in Sills & Bailyn (1999). 
Since the observed distributions of blue stragglers in the color-magnitude diagram are already 

extended to cooler temperatures than the predicted distributions, the increase in blue straggler 
main-sequence lifetime will exacerbate the lack of agreement between observations and theory. It 
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is possible that the width of the observed distribution is a selection effect, caused by the blended 
light of chance superpositions or true binary stars. More careful observations of blue straggler 
populations in globular clusters may clear up this particular mystery. 

While none of the stars investigated in this work arc fully mixed, rotation docs cause helium to 
be mixed to the surface. Rotating blue stragglers tend to be much brighter than their non-rotating 
counterparts, and much hotter. Therefore, rotation could be responsible for anomalously bright 
blue stragglers, such as the central 6 blue stragglers in NGC 6397 (Lauzeral et al. 1992) or the 
brightest three blue stragglers in M3 (Ferraro et al. 1997). Detailed comparisons between observed 
distributions and theoretical distributions based on rotational tracks are necessary to determine if 
there are two distinct populations of collisional blue stragglers in globular clusters, or if all blue 
stragglers can be explained by a population with a range of collisional impact parameters. 

The only blue straggler in a globular cluster to have its rotation rate measured is BSS 19 in 47 
Tucanae (Shara, Saffcr, & Livio 1997). This star is rotating at t;sini = 155 ± 55 km s^^, which is 
consistent with the rotation rates of the collision products presented in this work. The position of 
this star in the HR diagram is shown in Figure 7. The rotation velocities of the collision products 
in the vicinity of this star range from about 230 km to 400 km s"^, which is consistent with 
the BSS 19 data point and a reasonable inclination angle i between ~ 20° and ~ 40°. Observations 
of more rotation rates of blue stragglers in globular clusters are necessary to resolve some of the 
issues involved in studying off-axis collision products. However, we can say that the conclusion of 
Shara, Saffer, &; Livio (1997) that such a rapid rotation rate necessarily indicates a binary merger 
formation for BSS 19, rather than a collisional origin, is currently unwarranted. This conclusion was 
based on the assumption of Leonard Sz Livio (1995) that collisional blue stragglers would develop 
convective envelopes shortly after the collision, and would then lose angular momentum through a 
magnetic wind. As shown in Paper I and confirmed in this work, collision products do not develop 
significant surface convection zones, and therefore do not suffer from standard magnetic braking. 

The other blue stragglers with measured rotation rates arc found in the old open cluster M67 
(Peterson, Carney & Latham 1984). These stars are rotating more slowly than the blue straggler 
in 47 Tuc, and much more slowly than the models presented here. A large fraction of these blue 
stragglers are in binary systems. Could these two circumstances be related? Blue stragglers are 
thought to be formed through mass transfer in a binary system as well as through direct stellar 
collisions. Detailed calculations will be necessary to determine if binary merger products have 
different rotational properties than stellar collision products, and if we can distinguish between the 
two blue straggler creation mechanisms on the basis of their current rotation rates. 
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Fig. 1. — Profiles of the entropy A, pressure P, density p and radius r as a function of enclosed 
mass m for realistically modeled 0.6Mq (dashed curve) and O.SM© (dotted curve) parent stars. 
Units are cgs. 
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Fig. 2. — Snapshots of those SPH particles near the orbital plane (i.e., within two smoothing 
lengths of z = 0) for parabolic collision between parent stars of masses Mi = O.SMq 

and M2 = O.GMq at a pericenter separation Vp = 0.25(i?i + R2). Colors are used to indicate from 
which parent star the SPH particle originated, as well as its initial entropy A: green and red points 
are for particles originating in the 0.6Mq star, while yellow and blue indicate particles from the 
O.8M0; furthermore, green and yellow correspond to particles with an entropy A < 10^^ (cgs), 
while red and blue correspond to particles with larger entropies. (This figure can be downloaded 
from http://www.astronomy.ohio-state.edu/ asills/papers/rotBS.html) 
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Fig. 3. — Constant density, specific angular momentum and entropy contours at the end of each 
simulation in slices containing the remnants' rotation axes. Here rcyi is the cylindrical radius 
measured from the rotation axis, in units of the radius of a turnoff mass star. The closed loops that 
extend to both the right and left halves of each plot correspond to the isodensity surfaces enclosing 
15%, 30%, 45%, 60%, 75% and 90% of the remnant mass; the thick outermost bounding curve 
correspond to an enclosed mass fraction m/M = 0.95. The left half of each plot shows constant 
entropy contours, with a linear spacing of 1.04 x 10^^ (cgs), while the right half of each plot shows 
constant specific angular momentum contours, with a linear spacing of 2.66 x 10^^ cm^ s~^. 
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Fig. 4. — Angular velocity uj as a function of the enclosed mass fraction m/M in the final merger 
remnants. The points represent particle values from an SPH simulation, with only half of the 
particles displayed in cases e, f and k. The solid lines show the oo profiles of the initial YREC 
models, generated by a procedure that constrains the entropy and specific angular momentum 
profiles of each model to be that given by the SPH results. Units of u are rads~^. 
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Fig. 5. — Evolutionary tracks for case k with the initial angular velocity loq divided by 5 (dotted 
line), 10 (short dashed line), 100 (long dashed line) and 1000 (solid line line). The end of the 
collision is marked with a 'C, the start of central hydrogen burning with 'ZAMS', the turnoff with 
'TO', and the base of the giant branch with 'GB'. 
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Fig. 6. — Evolutionary track for case j' (solid line), and for the same case with the same initial 
structure profiles, but zero angular momentum (dotted line). The different evolutionary points are 
marked as in Figure 5. 
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Fig. 7. — Evolutionary tracks for case e (solid line), case f (dotted line), and case k (dashed line), 
with the initial rotational velocity divided by 5 in all cases. The different evolutionary points are 
marked as in Figure 5. The data point is BSS 19 in 47 Tuc from Shara, Saffer, & Livio (1997). 
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Fig. 8. — Surface helium abundance as a function of age for each of the four cases: cases c, f and 
k with the initial rotational velocity divided by 5 in each case, and case j'. The ages of the main 
sequence, the turnoff and the base of the giant branch are marked. 
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Fig. 9. — Surface rotational velocity as a function of time for four cases: case e with wq/S, case f 
with wo/S, case k with wq/S, and case j'. 
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Tablo 1. SUMMARY OF COLLISIONS 



Case 


Ml 


Ma 


rp 


N 


*/ 


J 


T/\W\ 


Mr 




Vy 








[Ri + R2\ 




[hours] 


[g cm^ 




[Mq] 


[km s"-*-] 


[km s-l] 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


(9) 


(10) 


(11) 


e 


0.8 


0.6 


0.25 


1.05 X 10^ 


11.10 


2.0 X 10^1 


0.101 


1.359 


-1.8 


-2.7 


f 


0.8 


0.6 


0.50 


1.05 X 10^ 


24.64 


2.8 X 10^1 


0.119 


1.380 


-0.4 


-1.2 


k 


0.6 


0.6 


0.25 


9 X 10" 


11.10 


1.3 X 10^1 


0.085 


1.158 


0.0 


0.0 


j' 


0.6 


0.6 


0.01 


1.6 X 10^ 


9.25 


2.1 X 10^° 


0.005 


1.147 


0.0 


0.0 



Table 2. EVOLUTIONARY RESULTS 



Case 


Vrot (TO) 


age at ZAMS 


age at TO 


age at GB 




(km s"-*-) 


(Gyr) 


(Gyr) 


(Gyr) 


e 


344 


0.012 


0.58 


0.87 


f 


372 


0.040 


2.79 




k 


353 


0.077 


4.19 


4.38 


j' 


272 


0.091 


3.66 


4.03 



